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SUMMARY 

Two  refinements  in  the  application  of  numerical  integration,  alternative  sampling 

points  and  reduced  integration,  liave  been  investigated  for  a Galerkin  finite  element 
formulation  of  a representative  How  problem.  By  sampling  the  solution  at  the  Gauss 
points  a significant  improvement  in  accuracy  is  achieved.  The  accuracy  gain  is  lost  if 
the  solution  at  the  Gauss  points  is  extrapolated  to  the  edge  of  the  element.  Consider- 
ation of  a one-dimensional  problem  suggests  that  the  use  of  reduced  integration  is 
equivalent  to  fitting  the  equation  residual  in  the  least-squares  sense  over  each  element. 
The  employment  of  reduced  integration,  rather  than  exact  integration,  for  incompress- 
ible, inviscid  flow  about  a two-dimensional  cylinder  has  produced  solutions  that  are  ten 
times  more  efficient  if  quadratic  rectangular  elements,  either  Serendipity  or  Lagrange, 
are  used.  The  utilization  of  linear  rectangular  elements  has  caused  a smaller  improve- 
ment. The  improvements  associated  with  the  introduction  of  reduced  integration  are 
independent  of  grid  refinement.  The  use  of  reduced  integration  and  triangular  ele- 
ments, with  both  linear  and  quadratic  shape  functions,  lias  produced  no  significant 
improvement.  The  results  herein  reinforce  the  previously  published  conclusion  that 
the  quadratic,  rectangular.  Serendipity  element  is  the  most  efficient  element  for  flow 
problems  of  this  type,  a 
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I.  INTRODUCTION 

The  application  of  the  finite  element  method  to  fluid  flow  problems  has  been  a fairly  recent  occurrence.  Most 
applications  so  far  have  been  aimed  at  demonstrating  the  possibility  of  obtaining  solutions  to  fluid  flow  problens  by 
the  finite  element  method  rather  than  obtaining  the  most  efficient  solution.  As  a consequence  a number  of  the  ^ 

refinements  that  have  been  employed  in  structural  applications  of  the  finite  element  method  to  improve  efficiency  have 
yet  to  be  tried  on  fluid  flow  applications.  Computational  efficiency  is  determined  by  the  execution  time  required  to 
achieve  a predetermined  accuracy. 

The  purpose  of  this  report  is  to  consider  two  such  refinements,  alternative  sampling  points  and  reduced  integration,  | 

applied  to  a representative  fluid  flow  problem.  Both  these  refinements  relate  to  the  numerical  integration  that  is 
necessary  to  convert  the  governing  ^tial  differential  equations  into  governing  algebraic  equations.  These  concepts  , 

have  been  explored  in  such  applications  as  the  analysis  of  plates  and  shells(ref.  1 and  2),  the  smoothing  of  discontinuous 
stresses(ref.3)  and  the  use  of  “thick  shell”  programs  for  solving  thin  shell  problems(ref.4);  a significantly  improved 
efficiency,  in  the  above  sense,  has  been  reported  in  many  cases. 

The  representative  fluid  flow  problem  is  the  flow  about  a two-dimensional  circular  cylinder  of  an  incompressible, 
inviscid  fluid.  This  problem  is  chosen  because  it  possesses  an  exact  solution  which  permits  a direct  measure  of  the 
accuracy.  A Galerkin  finite  element  formulation  in  primitive  variables  has  been  used  as  this  is  more  representative 
of  general  fluid  flow  problems  than  is  the  use  of  a variational  principle.  Previously  the  same  model  problem  and  i 

formulation  have  been  used  to  systematically  assess  the  relative  merits  of  various  elements  and  shape  functions(ref.S). 

The  plan  of  this  report  is  as  follows.  In  Section  2 some  background  developments  of  alternative  sampling  points  and  : 

reduced  integration  are  introduced.  In  Section  3 results  for  various  sampling  points  with  different  orders  of  numerical 
integration  are  presented  and  compared.  Additional  results  are  given  in  Section  3 which  illustrate  an  inaeased  accuracy  1 

for  the  same  execution  time  when  reduced  integration  is  used  instead  of  full  integration.  The  same  elements  and  shape  ' 

functions  that  were  compared  for  full  integration  in  reference  5 are  compared  for  reduced  integration  in  Section  3. 

This  comparison  considers  both  coarse  and  refined  grids.  In  Section  4 the  results  of  Section  3 are  interpreted  in  the 
light  of  the  background  material  given  in  Section  2.  Section  5 contains  an  explicit  list  of  the  conclusions  of  this  study. 


2.  BACKGROUND 

In  this  section  a little  of  the  background  development  of  alternative  sampling  points  and  reduced  integration  is  set  down. 

Although  both  concepts  have  been  established  qualitatively, published  material  indicates  that  very  few  quantitative 
results  are  applicable  to  other  than  specific  cases. 

2.1  Alternative  sampling  points 

A typical  finite  element  formulation  (e.g.  reference  S)  converts  a problem  of  continuous  variables  (e.g.  u,  v) 
governed  by  partial  differential  equations  into  a problem  of  discrete  variables  (e.g.  nodal  values  u governed  by 
algebra!  c equations.  This  is  done  by  arbitrarily  approximating  the  continuous  variables  by  low  order,  piecewise 
polynomials  in  which  the  nodal  values  are  the  unimown  coefficients.  To  form  the  algebraic  relations  typically 
a weighted  form  of  the  governing  equations(ref.6),  with  the  polynomial  approximations  substituted,  is  integrated 
over  the  whole  domain.  If  the  integration  is  performed  analytically  the  only  approximation  present  is  associated 
with  forcing  the  solution  to  behave  locally  as  though  it  were  a low-order  polynomi  al,  assuming  an  isoparametric 
formulation  is  not  used. 

If  the  integration  is  performed  numerically  additional  errors  may  be  introduced.  Numerical  integration  implicity 
fits  a low  order  function,  typically  a polynomial,  to  the  data  and  performs  the  integration  analytically.  The  low 
order  function  fits  the  data  by  requiring  an  exact  match  at  certain  predetermined  points.  The  numerical  inte-  i 

gration  is  then  a weighted  sununation  of  the  integrands  evaluated  at  the  predetermined  points.  Normally  Gauss  ’ 

quadrature  is  used  because  with  n points  it  is  capable  of  integrating  exactly  integrands  of  polynomial  form  up  to  ; 

order  2n  - l(ref.7). 

Thus  normally  the  order  of  the  numerical  integration  is  chosen  to  be  sufficiently  high  that  the  integration  is 
exact  for  the  order  of  the  finite  element  representation  used.  ! I 

For  problems  that  can  be  solved  by  formulating  a variational  principle  it  has  been  shown(ref.4)  that  at  the  j 

Gauss  points  the  differential  expressions  in  the  integrand  are  closer  to  the  exact  values  then  elsewhere.  For  1 

example,  for  a cantilever  beam  subject  to  a distributed  load  the  stresses  are  determined  very  accurately  if  sampled  I 

at  the  Gauss  points  but  are  very  inaccurate  and  discontinuous  if  sampled  at  the  nodal  points.  i 

In  the  present  problem  alternative  points  within  the  element  are  considered  to  see  if  consistently  better 
agreement  with  the  exact  solution  can  be  obtained.  . | 
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2.2  Reduced  integration 

Reducing  the  order  of  integration  below  that  required  to  make  the  integration  exact  necessarily  introduces  some 
additional  error.  Cleaily  if  the  order  is  reduced  too  much  convergence  may  no  longer  be  achieved.  The  minimum 
order  of  integration  for  convergence  is  that  which  permits  exact  integration  of  the  element  volume  (or  area  in  two 
dimensions)  when  an  isoparametric  formulation  is  used(ref.8). 

Ilie  optimum  order  of  integration  depends  on  the  order  of  the  finite  element  representation,  p,  and  the  order  of 
differentiation,  d,  in  the  governing  equation.  For  problems  of  a variational  nature  such  as  the  plane  strain  problem 
it  has  been  deduced(ref.4)  that  for  p = 1 and  d = 1 , 1 point  integration  formulae  should  be  optional  while  for 
p = 2 and  d = 1 , 4 point  integration  formulae  would  be  required  if  an  element  is  considered  in  isolation. 

If  the  stiffness  matrix  is  considered  as  a whole,  reduced  integration  may  be  viewed  from  a different  perspective. 
Since  the  Galerkin  method  is  an  example  of  the  Method  of  Weighted  Residuals  it  is  instructive  to  derive  the 
governing  algebraic  equations  from 


!l 


N. . R . dx  dy  = 0, 


(1) 


where  R is  the  residual  of  one  of  the  equations  after  the  analytic  representation  in  terms  of  the  nodal  unknowns 
has  been  substituted  (see  reference  6 for  more  details).  N.  is  the  shape  function  appropriate  to  the  i^^  non,.. 

If  the  integration  in  equation  ( I ) is  carried  out  numerically  then  the  result  is 


Wk.N.(x^.yk).R(x^,yk)  = 0. 


(2) 


is  the  weight  attached  to  the  k function  evaluation  point.  M is  the  number  of  function  evaluation  points, 

i.c.  the  produce  of  m function  evaluation  points  per  clement  and  c contributing  elements.  Because  of  the  nature 
ol  N.  only  the  four  rectangular  elements  surrounding  the  i*^  node  will  contribute  to  the  summation.  For  example 

al  node  I in  figure  I , only  elements  w,  x,  y and  /.  contribute.  Thus  in  any  one  equation  M <4  m. 

The  residual  R could  depend  on  any  of  the  n nodal  unknowns,  q.,  in  the  domain.  Thus 


Typically 


II 

3N. 


(3) 


j'^'k-^k' 


For  convenience,  elements  are  assumed  to  be  perfectly  regular.  The  symbol  n is  the  product  of  the  number  of 
nodes  n^  and  the  number  of  unknowns  per  node  which  in  this  case  is  2.  Certain  of  the  q.  will  be  known  and 

hence  equation  (2)  with  (3)  substituted  can  be  written. 


M L_ 

i^i  >1 

where  B.  represents  the  weighted  summation  associated  with  the  known  q..  Th®  problem  now  has  L unknowns 
and  L equations  like  (4)  must  be  written  down  in  order  to  obtain  a unique  solution  for  the  q.’s.  It  is  also 

apparent  that  for  a solution  to  t>e  possible,  the  total  number  of  function  evaluations  M.^  must  not  be  less  than  the 
number  of  unknowns  L 


-3- 
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Using  exact  numerical  integration,  is  much  larger  than  L,  using  reduced  integration  is  much  closer  to  L 

although  still  greater.  There  is  a possibility  that  if  L is  too  close  to  a subsystem  of  [q]  might  be  completely 

determined  from  a few  elements;  this  is  clearly  a singular  situation. 

If  the  order  of  summation  in  equation  (4)  is  reversed  the  familiar  stiffness  equation  is  obtained, 


L 

E 

It 

(5) 

M 

where  K.. 

U 

II 

N.(x^,y^.Ri.(x^,y^). 

Consequently  for  a solution  using  reduced  integration  to  lead  a non-singular  stiffness  matrix  [K] , it  is 
necessary  that 


I 

i 


e . m > L (6) 

I j 

The  use  of  reduced  integration  implies  a fitting  of  the  integrand  with  a lower  order  polynomial  than  that  I ‘ 

implied  by  the  order  of  the  finite  element  approximation.  Hinton  and  Campbell(ref.3)  prove  that  for  two- 

* dimensional  parabolic  isoparametric  rectangular  elements  applied  to  a typical  variational  formulation  of  the 
finite  element  method,  the  use  of  reduced  integration  is  equivalent  to  performing  the  integration  exactly 
after  applying  a least  squares  bilinear  fit  separately  to  each  term  of  the  integrand. 

• Thus  since  the  finite  element  formulation  constrains  the  solution  to  locally  follow  a low-degree  polynomial  the 
use  of  reduced  integration  may  be  thought  of  as  relaxing  some  of  the  constraints  on  the  solution  at  the  expense 
of  introducing  additior:al  errors  by  not  carrying  out  the  integration  exactly.  The  empirical  evidence  clearly 
indicates  that  the  g^ins  associated  with  relaxing  the  constraints  far  outweigh  the  losses  associated  with  executing 
the  numerical  integration  less  accurately. 

A fuller  discussion  of  specific  applications  of  reduced  integration  may  be  found  in  reference  4 and  a discussion 
of  reduced  integration  in  a wider  context  may  be  found  in  reference  9. 

2.3  A worked  example 

A simple  example  in  one  dimension  will  be  considered  in  order  to  illustrate  how  and  why  reduced  integration 
works. 

The  governing  equation  is 


; Uy)  = ^-y  = o,  (7) 

1 

subject  to  the  boundary  condition  y = 1 at  x = 0.  It  is  of  historical  interest  that  this  problem  was  used  by 
. Duncan(ref.lO)  to  illustrate  the  classical  Galerkin  method  (see  also  reference  6).  This  problem  has  the  exact  solution 

I X 

r y = e - 

. A solution  is  sought  for  (Xx  < I and  five  nodal  values  y , at  equal  intervals  of  x,  are  introduced.  The  domain 

1 is  split  into  two  elements,  (Xx  <0.5  and  O.S  <x  < 1 , and  quadratic  Lagrange  shape  functions  are  used  to 

I determine  the  local  variation  of  y. 


From  the  boundary  condition  7i  = 1 ; hence  only  four  independent  equations  are  required  to  uniquely 
determine  the  unknown  y.’s.  The  analytic  representation  of  y is  given  by 
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Substituting  this  into  equation  (7)  produces  a residual,  R(x),  given  by 

5 


Forming  the  integral  of  the  weighted  residual  over  the  domain  (see  reference  6)  gives 

1 

[ W.  R(x)  dx  = 0. 


For  the  Galerkin  method  W.  = N.  and  equation  (10)  becomes 


N.  I - N.  1 dx.  y.  = 

‘ I J J 


Setting  i = 2 to  S gives  sulTicicnt  algebraic  equations  of  the  form 

5 


a.. . y.  =0. 
ij  J 


to  solve  for  the  unknowns  : y2,  ya,  ya  and  ys- 


The  coefficients  a„  are  given  by 


_j  - N dx 
dx  ^ 


and  are  to  be  evaluated  numerically.  The  integral  is  evaluated  for  each  element  separately.  If  a three  point  Gauss 
quadrature  formula  is  used  the  integrations  are  performed  exactly  and  the  corresponding  solution  for  the  nodal 

values  y.  arc  shown  in  Table  I . 

J 


tabu:  I.  COMPARISON  OF  SOLUTIONS  FOR  ^ - y = 0 

dx 


Approximate  solution  (nodal),  y 


Exact 

numerical 

integration 

9" 

Reduced 

numerical 

integration 

y" 

Local  L.S. 
fit  of 
residual 

y'* 

Traditional 
Galerkin 
(quart  ic 
polynomial) 

1.0000 

1.0000 

1.0000 

1.0000 

1.2707 

1 .2838 

1.2838 

1.2840 

1.6403 

1.6486 

1.6486 

1.6488 

2.0990 

2.1165 

2.1165 

2.1170 

2.6938 

2.7180 

2.7180 

2.7182 

Exact 

solution 

X 

e 


1.0000 

1.2840 

1.6487 

2.1170 

2.7183 
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If  a two  point  Gauss  quadrature  (reduced  integration)  formula  is  used  the  integrations  are  not  formed  exactly; 
the  corresponding  results  for  the  nodal  values  are  shown  in  Table  I . It  is  self-evident  that  the  use  of  reduced 
integration,  rather  than  exact  integration,  has  produced  a very  much  more  accurate  solution.  It  would  be  possible 
to  construct  a local  least -squares  fit  to  the  term 


in  equation  (1 3)  and  carry  out  the  integration  exactly.  If  this  were  done  the  results  would  be  identical  to  those 
shown  in  Table  1 under  reduced  integration. 

This  is  because  the  least -squ^re^  fit  of  order  n - 1 always  cuts  the  original  curve  of  order  n at  the  n Gauss  points. 
Thus  the  least -squares  fit  of 


‘JN.  M 
dx 


over  each  element  would  be  a straight  line  cutting  the  original  quadratic  function  at  the  two  Gauss  points. 
Since  the  weight  function  N.  is  quadratic,  the  total  integration,  over  each  element,  of 


! K !™|  .N. 

Jo  ' Isf  ' 


could  be  carried  out  exactly  using  a two  point  Gauss  quadrature  formula.  Thus  the  sampling  values  of  the  integrand 
would  be  exactly  the  same  as  for  the  reduced  integration  of  the  original  expression  in  equation  (13). 

Also  shown  in  Table  1 are  results  for  a traditional  Galerkin  formulation(ref.6)  in  four  unknowns.  It  is  generally 
found  that  because  the  analytic  representation  and  the  weighting  functions  are  allowed  to  span  the  complete  domain 
superior  results  are  obtained(ref.6)  with  the  traditional  Galerkin  method  than  with  a fmite  element  or  finite 
difference  method  with  the  same  number  of  unknowns.  Table  1 indicates  that  the  use  of  reduced  integration  has 
produced  results  that  arc  almost  as  good  as  using  a traditional  Galerkin  formulation. 

The  question  remains:  why  does  the  use  of  reduced  integration  give  more  accurate  results  than  the  use  of  exact 
integration?  The  answer  can  be  found  by  pursuing  the  equivalence  between  the  use  of  reduced  integration  and 
'he  least  squares  fitting  of  the  expression 


‘IN.  M 

-'•"j 


in  equation  (13)  and  subsequent  integration.  If  equations  (13),  (1 2),  (1 1)  and  (10)  are  considered  it  is  clear  that 
the  least-squares  fit  of  the  expression 

15  ",] 

^ dx 

over  each  element  is  equivalent  to  a least -squares  fit  over  each  element  of  the  residual  in  equation  (10). 

The  basis  for  the  weighted  residual  method(ref.6)  is  that,  since  the  introduction  of  an  analytic  representation 
produces  a non-zero  residual  equation  throughout  the  domain,  a reasonable  approximate  solution  will  be  obtained 
if  the  residual  is  made  to  equal  zero  in  some  global  sense  e.g.  equation  (10).  Ideally  the  residual  should  be  zero 
everywhere.  If  the  weight  function  is  predominantly  of  one  sign  the  variation  of  R may  be  expected  to  include 
some  changes  of  sign  if  equation  (10)  is  to  be  satisfied.  Therefore  a low  order  least-squares  fit  of  R whose 
weighted  integral  over  the  domain  is  zero  is  likely  to  have  local  magnitudes  of  R^^  that  are  less  than  the  local 

magnitudes  of  R,  and  hence  produce  a solution  closer  to  the  exact. 
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This  is  illustrated  in  figure  2.  From  the  solution  using  exact  numerical  integration  the  expression 


dN.  . , ‘|  . ei 

J • N.  . y. 

.*r  'I 


has  been  plotted.  R^.  is  the  equation  residual  corresponding  to  the  approximate  solution,  j^y.  J R . is 
quadratic  in  x and  discontinuous  at  the  element  boundary  (x  = 0.5).  If  the  expression 


in  equation  (1 3)  is  fitted  by  a linear  least-squares  curve  a pseudo-residual,  based  on  y^  * , can  be  constructed  as 


j=l  ^ 

Thb  b also  plotted  in  figure  2.  It  is  apparent  that  the  effect  of  the  least -squares  fit  has  not  greatly  reduced  the 
size  of  the  residual.  However  the  residual  based  on 


dx  Jls 


“Is  “Is 

and  the  corresponding  solution  y^  , is  nmch  smallor.  y/  is  obtained  by  solving 


(a^.li^.y.  = 0, 


Rj^  is  given  by 


Rj^  is  actually  zero  but  has  been  plotted  in  figure  2 in  order  that  it  will  show  up.  Because  of  the  local  least-squares 
fit  R|j  is  a linear  function  of  x.  R|^  is  identical  with  the  residual  obtained  from  a reduced  integration  formulation,  R 
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The  lower  ihc  order  of  the  least-squares  Tit  the  closer  the  residuals  should  be  to  zero.  However  as  noted 
previously  the  total  number  of  function  evaluations  in  the  domain  must  not  be  less  than  the  number  of  unknowns. 
In  the  present  example  the  number  of  unknowns  is  four,  there  are  two  elements,  therefore  a two  point  Gauss 
quadrature  formula  is  the  lowest  order  possible. 

This  interpretation  of  reduced  integration  clearly  points  towards  the  use  of  a least -squares  approximation  to  the 
method  of  weighted  residuals  or  an  MWLSR  as  a valuable  technique  in  its  own  right.  The  “magic”  in  the  use  of  a 
lower  order  Gauss  quadrature  formula  is  that  it  permits  a shortcut  to  be  taken  so  that  the  least-squares  approxim- 
ation to  the  residual  is  set  up  implicitly.  Thus  the  use  of  Gauss  quadrature  formulae  is  convenient  but  not  essential 
to  the  success  of  a solution  after  approximating  the  residual  in  the  least -squares  sense. 


3.  RESULTS 


The  model  fluid  flow  problem  considered  in  this  report  is  inviscid,  incompressible  flow  about  a two-dimensional 
circular  cylinder;  the  flow-field  is  represented  schematically  in  figure  1 . The  governing  equations,  for  inviscid, 
incompressible  flow  in  two  dimensions,  are  taken  to  be 


^ -r  ^ = 0 
3x  dy 


(19) 


and 


dy  ax 


(20) 


After  application  of  a Galerkin  finite  element  formulation  the  governing  partial  differential  equations  can  be  reduced 
to  algebraic  equations  of  the  form 


j=l 

n 


a . . u.  -r  \ b.. . V.  = 0,  j = 1 , n 
'J  J 'J  J ^ 

j=l 

n 


(21) 


b.. . u.  - ) a.. . V.  = 0,  j = I,  n, 

ij  J / >J  J 


(22) 


j=l 


j=l 


where 


'dN. 


dN. 


a„  _j  . N.  dxdyandb.j  “yy  _J  - N,  dxdy. 


The  symbols  u.,  v.  represent  the  nodal  values  of  u and  v,  n is  the  total  number  of  nodal  unknowns  and  N.  is  the  shape 
11 

function  associated  with  the  i'"  node.  Further  details  may  be  found  in  reference  5. 


A computational  solution  is  sought  within  the  region  ABCD  shown  in  figure  1 . The  nodal  points  and  elements  are 
defined  on  a polar  grid  and  an  isoparametric  formulation  is  used  to  connect  this  to  a cartesian  grid. 

All  results  presented  in  this  section  are  for  the  variation  of  the  tangential  velocity  component  at  the  body  surface  with 
angular  position.  As  in  reference  S it  is  felt  that  this  represents  a stringent  test  of  the  computational  results. 

In  order  to  have  a single  number  represent  the  accuracy  of  the  computational  results,  the  root  mean  square  difference, 
a,  between  the  finite  element  solution  and  the  exact  solution  is  defined  as  follows 
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(23) 


where  q.^,  is  the  finite  element  solution  for  the  tangential  velocity  component.  The  symbol,  is  the  exact  solution 

for  the  tangential  velocity  component  and  N is  the  number  of  nodes  between  6 = 0°  and  90°  (B  and  C in  figure  I). 
The  symbol,  o may  be  thought  of  as  the  average  difference  between  the  computational  and  exact  solutions. 


3.1  Alternative  sampling  points 

From  Section  2.1  it  is  to  be  expected  that  the  solution  at  the  Gauss  points  will  be  closer  to  the  exact  solution 
calculated  at  the  Gauss  points  than  the  nodal  solution  is  to  the  exact  solution  at  the  nodes.  For  practical  purposes 
the  solution  at  the  body  surface  is  of  interest.  The  body  surface  coincides  with  the  boundary  of  a number  of 
adjacent  elements.  Thus  nodal  values  at  the  body  surface  are  available  buthe  Gauss  points  do  not  coincide  with 
the  body  surface. 

Using  the  analytic  representation  implicit  in  the  flnite  element  formulation  the  solution  at  the  Gauss  points  can 
be  obtained.  Using  just  the  solution  at  the  Gauss  points  an  extrapolation  is  made  onto  the  body  surface  as  shown 
in  figure  3.  This  gives  the  solution  C.  The  nodal  solution  B is  interpolated  along  the  body  surface  to  give  the 
solution  D which  is  at  the  same  point  as  C.  This  process  is  repeated  for  2 x 2 and  3x3  integration  formulae. 


Ii 

i' 
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TABLE  2.  R.M.S.  ERRORS  FOR  ALTERNATIVE  SAMPUNG  POINTS  AND  DIFFERENT 
NUMERICAL  INTEGRATION  FORMULAE 


Solution  A 

Solution  B 

Solution  C 

Solution  D 

Integration 

formula 

R.M.S.  difference 
at  Gauss  points 

R.M.S.  difference 
at  nodal  points 

Solution  A 
extrapolated  to 
the  surface 

Solution  B 
interpolated  to 
same  position  as 
Solution  C 

4 point 

0.085 

0.017 

0.016 

0.015 

9 point 

0.009 

0.062 

0.066 

0.058 

The  results  for  the  r.m.s.  differences  are  shown  in  Table  2.  These  results  were  obtained  using  second  order 
rectangular  elements  of  the  Serendipity  type  with  149  unknowns  in  the  whole  flow-field.  The  actual  differences 
between  the  computed  and  exact  solutions  are  plotted  against  angular  position  in  figures  4 and  S corresponding 
to  the  different  integration  formulae.  On  each  of  figures  4 and  S the  solutions  B,  C and  D (figure  3)  are  shown. 

Examination  of  Table  2 indicates  that  the  solution  at  the  Gauss  points  is  better  for  the  3 x 3 integration  formulae, 
i.e.  the  lowest  order  that  still  performs  the  numerical  integration  exactly.  This  solution,  using  a 3 x 3 integration 
formula  is  superior  to  the  nodal  point  solution  using  reduced  integration.  Once  the  Gauss  point  solution,  for  the 
3x3  integration  formula,  is  extrapolated  to  the  surface  (Solution  C in  Table  2)  the  increas^  accuracy  of  the 
Gauss  point  solution  is  lost.  In  fact  all  solutions  at  the  body  surface  are  of  comparable  accuracy. 

Tlie  use  of  a reduced  integration  formula  (2  x 2)  produces  a relatively  inaccurate  solution  at  the  Gauss  points 
(Solution  A)  but  a very  accurate  solution  when  extrapolated  to  the  body  surface  (Solution  C).  The  interpolation 
of  the  nodal  point  solution  (Solution  D)  to  the  same  position  as  the  extrapolated  Gauss  point  solution  (Solution  C) 
has  always  produced  a slightly  better  solution. 

Both  integrating  formulae  (Hgures  4 and  S)  produce  solutions  which  show  sonae  disagreement  between  solutions 
C and  D at  the  forward  stagnation  point  but  good  agreement  as  the  shoulder  point  (6  = 90°)  is  approached. 

The  general  conclusion  of  these  results  is  that  the  use  of  alternative  sampling  points  (Solutions  C and  D)  produces 
only  marginal  improvement  in  accuracy  when  compared  with  the  large  improvement  in  accuracy  at  the  nodes 
associated  with  the  use  of  reduced  integration.  Comparable  results  were  obtained  when  a refined  grid  was 
considered. 
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3.2  Reduced  integration  vs.  conventional  integration 

Solutions  have  been  obtained  using  a second  order  rectangular  element  of  the  Serendipity  type  with  a 2 x 2 
integration  formula  (reduced  integration)  and  a 3 x 3 integration  formula  (full  integration).  The  results  for  a 
coarse  grid  (149  nodal  unknowns  in  the  flow  field)  are  shown  in  figure  6.  Because  of  the  use  of  a dummy  element 
to  carry  out  the  numerical  integration  (see  reference  5)  the  use  of  reduced  integration  has  not  produced  any 
significant  improvement  in  execution  time  (as  suggested  in  reference  4).  Therefore  any  execution  time  shown  is 
nominally  the  same  irrespective  of  the  order  of  numerical  integration. 

The  results  for  the  coarse  grid  indicate  that  the  solution  using  reduced  integration  is  more  than  three  times  as 
accurate  as  the  solution  using  full  integration.  An  examination  of  figure  6 indicates  that  the  full  intepation 
solution  always  lies  below  the  exact  solution  whereas  the  reduced  integration  solution  straddles  the  exact  solution. 
Thus  an  additional  advantage  of  the  reduced  integration  solution,  not  apparent  from  the  r.m.s.  error  results,  is 
the  possibility  of  fitting  a low  order  least  squares  curve  through  the  reduced  integration  results  to  obtain  even 
better  agreement  with  the  exact  results. 

The  possibility  exists  that  the  improvement  produced  by  using  reduced  integration  will  lessen  as  the  grid  is 
refined.  Results  for  a moderately  refined  grid  (299  nodal  unknowns  in  the  flow  field)  are  shown  in  figure  7. 

As  is  apparent  the  reduced  integration  solution  is  still  better  than  three  times  as  accurate  as  the  full  integration 
solution.  A further  refinement  in  the  grid  (with  582  nodal  unknowns  in  the  flow  field)  also  produced  a reduced 
integration  solution  that  was  better  than  three  times  as  accurate  as  the  full  integration  solution.  Because  both 
solutions  are  very  close  to  the  exact  solution  this  case  has  not  been  plotted. 

A better  appreciation  of  the  improvements  associated  with  the  use  of  reduced  integration  may  be  obtained  by 
considering  figure  8.  The  reduced  integration  results  correspond  to  those  plotted  in  figure  6;  the  full  integration 
results  have  been  obtained  by  refining  the  grid  until  a solution  of  the  same  total  accuracy  is  produced.  It  may  be 
noted  that  this  has  required  more  than  ten  times  the  execution  time.  By  plotting  the  tangential  velocity 
differences  against  angular  position  the  tendency  for  the  full  integration  solution  to  underestimate  the  exact 
solution  and  the  reduced  integration  solution  to  straddle  the  exact  solution  is  made  more  obvious. 

3.3  Reduced  integration  on  a coarse  grid 

fn  this  section  various  element/shape  function  combinations  are  compared  for  solutions  on  a coarse  grid  using 
reduced  integration.  The  same  element/sliape  function  combinations  are  considered  here  as  were  considered  for 
full  integration  in  reference  S,  namely; 

(1 ) linear  sliape  functioii  in  triangular  isoparametric  element 

(2)  quadratic  shape  function  in  triangular  isoparametric  element 

(3)  linear  shape  function  in  rectangular  isoparametric  element 

(4)  quadratic  shape  function  (Lagrange)  in  rectangular  isoparametric  element 

(5)  quadratic  shape  function  (Serendipity)  in  rectangular  isoparametric  element 

All  these  shape  functions  are  described  in  reference  1 1 . The  various  numerical  integration  formulae  used  with 
various  element/shape  function  combinations  are  described  in  Appendix  I.  However  the  results  presented  in 
this  section  cannot  be  compared  with  those  given  in  reference  S since  the  grid  used  here  is  coarser;  the  grids 
considered  in  reference  5 correspond  to  the  grids  used  in  Section  3.4. 

Results  for  linear  and  quadratic  shape  functions,  for  rectangular  elements,  are  compared  in  figure  9.  The 
numerical  integrations  for  the  quadratic  shape  functions  have  been  evaluated  with  a 4 point  formula  based  on 
evaluating  the  integrands  at  the  Gauss  points  (see  Appendix  I).  The  numerical  integrations  for  the  linear  shape 
functions  have  been  evaluated  with  a 1 point  formula  based  on  evaluating  the  integrands  at  the  mid  point  of  the 
elements  (see  Appendix  I).  Although  both  solutions  straddle  the  exact  solution  the  solution  using  a quadratic 
shape  function  is  noticeably  smoother  particularly  close  to  d = 90°.  As  is  apparent  from  Table  4 the  solution, 
using  quadratic  shape  functions  of  the  Serendipity  type,  is  approximately  twice  as  accurate  as  the  solution  using 
linear  shape  functions.  The  solution  using  quadratic  shape  functions  is  abo  more  economical,  mainly  because 
it  has  less  nodal  unknowns.  A comparison  of  Tables  3 and  4 indicates  that  the  use  of  reduced  integration  in  place 
of  exact  integration  (3x3  formula)  causes  a greater  improvement  in  accuracy  for  quadratic  shape  functions  than 
it  does  for  linear  shape  functions. 

Results  for  linear  and  quadratic  shape  functions,  for  triangular  elements,  are  shown  in  figure  10.  The  numerical 
integrations  for  the  quadratic  shape  functions  have  been  evaluated  using  a 4 point  formula  given  by  Cowper(ref.1 2). 
The  numerical  integration  for  the  linear  shape  functions  have  been  evaluated  with  a 1 point  formub  given  by 
Zienkiewicz(ref.I  I).  Examination  of  figure  10  and  Table  4 indicates  that  the  results  for  linear  elements  are 
closer  to  the  exact  solution  on  average  but  become  rather  poor  close  to  0 <=  90°.  The  solution  using  linear  shape 
functions  b exactly  the  same  whether  a I point  (Table  4)  or  7 point  (Table  3)  integration  formula  is  used. 
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This  indicates  that  use  of  a I point  integration  formula  is  exact  for  linear  triangular  elements  applied  to  the 
present  problem.  The  use  of  reduced  integration  for  a quadratic  triangular  element  produces  some  improvement 
in  accuracy  (Tables  3 and  4). 

Results  for  triangular  and  rectangular  elements,  with  linear  shape  functions  are  shown  in  figure  1 1 . Since  the 
results  with  triangular  elements  are  exact  with  the  use  of  a one  point  integration  formula  the  results  generally  lie 
beneath  the  exact  results  whereas  the  results  for  the  use  of  rectangular  elements  straddle  the  exact  solution.  Both 
sets  of  results  have  the  same  number  of  nodal  unknowns  and  hence  approximately  the  same  CPU  time  but  the 
solutions  using  rectangular  elements  are  approximately  twice  as  accurate  (see  Table  4). 

Results  for  triangular  and  rectangular  elements,  with  quadratic  shape  functions  are  shown  in  figure  1 2.  The 
numerical  integrations  for  the  quadratic  shape  functions  have  been  evaluated  with  a 4 point  formula  (see  Appendix  I). 
The  results  for  the  rectangular  elements  (Lagrange  shape  function)  lie  much  closer  to  the  exact  solution  (figure  12) 
particularly  close  tod  = 90°.  Both  solutions  required  approximately  the  same  CPU  time  but  the  results  using 
rectangular  elements  are  almost  three  times  as  accurate  on  average.  The  use  of  reduced  integration  for  quadratic 
rectangular  elements  of  the  Lagrange  family  has  produced  a two-fold  improveriKnt  in  accuracy  compared  with 
the  use  of  exact  integration  (Tables  3 and  4). 

Quadratic  rectangular  elements  of  the  Lagrange  and  Serendipity  families  have  been  compared  and  the  results 
are  shown  in  figure  1 3.  The  grids  and  number  of  unknowns  have  been  adjusted  to  give  approximately  the  same 
CPU  time.  As  can  be  seen  from  figure  13  the  resuhs  using  quadratic  shape  functions  of  the  Serendipity  family 
lie  substantially  closer  to  the  exact  solution  than  the  results  using  a Lagrange  shape  function.  The  average 
accuracy  using  Serendipity  shape  functions  is  approximately  three  times  greater  than  the  accuracy  using  Lagrange 
shape  functions  for  comparable  CPU  times  (Table  4).  Both  types  of  shape  function,  with  reduced  integration, 
produce  results  that  straddle  the  exact  solution  and  hence  could  be  made  to  produce  better  agreement  by  a low 
order  least  squares  fit. 

3.4  Reduced  integration  on  a moderate  grid 

The  same  element/shape  function  combinations,  as  set  out  in  Section  3.3,  are  here  compared  for  solutions  on 
moderate  grids  using  reduced  integration.  The  detailed  results  are  shown  in  figures  14  to  18  and  the  gross 
properties  of  the  solutions  are  listed  in  Table  6.  Comparative  results  using  exact  integration  are  given  in  Table  S. 
Many  of  these  results  have  been  obtained  using  the  same  grids  as  in  reference  5.  However  the  results  shown  in 
Table  S have  required  considerably  less  CPU  time  than  the  results  presented  in  reference  5.  This  improvement  in 
efficiency  is  due  to  system  changes  to  the  computer  (IBM  370/168)  used  to  obtain  the  results  rather  than  to 
improvements  in  this  program. 

Results  for  linear  and  quadratic  shape  functions,  for  rectangular  elements,  are  shown  in  figure  14.  The  solution 
using  quadratic  shape  functions  of  the  Serendipity  type  is  closer  to  the  exact  solution  particularly  close  to  0 ==  90°. 
However  the  general  character  of  both  solutions  is  similar  to  that  of  the  solutions  using  coarse  grids  shown  in 
figure  9.  Comparison  of  Tables  4 and  6 indicates  that  the  solution  using  quadratic  shape  functions  is  approxima- 
tely twice  as  accurate  as  the  corresponding  solution  on  a coarse  grid.  The  results  using  linear  shape  functions  are 
also  more  accurate  than  the  conesponding  results  on  a coarse  grid.  Also  both  linear  and  quadratic  solutions  are 
considerably  more  accurate  than  corresponding  solutions  (Table  S)  obtained  using  an  exact  integration  formula 
(see  Appendix  I)  on  a moderate  grid,  although  the  improvement  is  greater  for  quadratic  shape  functions. 

Linear  and  quadratic  shape  functions  used  with  triangular  elements  have  been  compared  and  the  results  are 
presented  in  figure  1 5.  The  results  are  of  approximately  the  same  accuracy  although  the  solution  using  quadratic 
elements  is  more  accurate  close  toff  90°.  Both  solutions  are  considerably  more  accurate  than  those  obtained 
using  a coarse  grid  (figure  1 0).  As  with  the  coarse  grid  solutions,  the  results  using  a linear  shape  function  are 
exactly  the  same  whether  a I point  or  7 point  integration  formula  is  used.  Comparison  of  Tables  5 and  6 
indicates  that  the  solution  using  quadratic  shape  functions  is  less  accurate  with  reduced  integration  than  with 
exact  integration. 

The  results  for  triangular  and  rectangular  elements  using  linear  shape  functions  are  shown  in  figure  16  and  it  is 
apparent  that  the  rectangular  element  solution  is  considerably  more  accurate.  Since  the  triangular  element 
solution  is  the  same  whether  the  integration  is  exact  or  reduced  and  since  both  elements  produce  results  of  equal 
average  accuracy  when  the  integration  is  exact,  any  superiority  of  the  rectangular  element  must  be  due  to  the 
use  of  reduced  integration. 

A comparison  of  results  for  the  use  of  triangular  and  rectangular  elements  with  quadratic  shape  functions  on  a 
moderate  grid  indicates  (figure  1 7 and  Table  6)  that  the  use  of  rectangular  elements  of  the  Lagrange  type  produces 
results  that  are  three  times  as  accurate  as  those  using  triangular  elements.  This  improvement  in  accuracy  is  of  the 
same  order  as  was  obtained  with  coarse  grids  using  both  exact  and  reduced  integration. 
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However  the  use  of  exact  integration  on  a moderate  grid  produces  more  accurate  results  with  a quadratic 
triangular  element  than  with  a quadratic  rectangular  element  of  the  Lagrange  family  (Table  S). 

Figure  1 8 shows  results  of  a comparison  of  quadratic  rectangular  elements  of  the  Serendipity  and  Lagrange 
families.  Both  sets  of  results  lie  close  to  the  exact  solution;  an  examination  of  Table  6 indicates  that  the 
solution  using  Serendipity  elements  is  more  accurate.  Both  solutions  are  considerably  more  accurate  than 
conesponding  solutions  using  exact  integration  (Table  S)  and  more  accurate  than  corresponding  solutions 
obtained  on  a coarse  grid  (figure  13  and  Table  4). 

3.5  Effect  of  grid  refinement 

Solutions  have  been  obtained  for  various  meshes  using  isoparametric  rectangular  elements  with  quadratic  shape 
functions  of  the  Serendipity  family  and  using  reduced  integration.  The  results  for  three  representative  cases 
are  shown  in  figure  19.  The  results  using  a coarse  mesh  required  the  solution  of  149  nodal  unknowns  and 
required  4.5  s CPU  time;  the  results  for  a moderate  mesh  required  the  solution  of  299  nodal  unknowns  and 
required  1 5 s CPU  time  and  the  results  for  a refined  mesh  required  the  solution  of  582  nodal  unknowns  and 
required  69  s CPU  time.  As  is  apparent  from  figure  19  all  three  solutions  straddle  the  exact  solution  and  so 
could  be  made  to  agree  with  the  exact  solution  even  better  by  curve  fitting  the  data  in  the  least-squares  sense 
with  a low  order  function.  These  results  indicate  that  the  accuracy  increases  with  the  square  root  of  the  CPU 
time  i.e.  to  double  the  accuracy  requires  four  times  the  CPU  time.  Examination  of  figure  19  indicates  that 
the  use  of  the  coarsest  grid  shown  produces  average  differences  that  are  less  than  1%  of  the  maximum  velocity 
whilst  use  of  the  most  refined  grid  produces  average  differences  that  are  less  than  0.2%  of  the  maximum  velocity. 


4.  DISCUSSION 

Since  flow  problems  normally  require  the  solution  on  the  element  boundary  it  appears  that  alternative  sampling  points 
within  the  element,  e.g.  the  Gauss  points,  offer  no  advantage  over  the  nodal  points.  As  is  apparent  from  Section  3.1 , 
if  the  solution  at  the  Gauss  points  is  extrapolated  to  the  element  boundary  the  accuracy  is  not  significantly  better  than 
obtaitied  at  the  nodes. 

This  section  will  be  mainly  devoted  to  a consideration  of  the  results  of  Sections  3.2  to  3.4  in  the  light  of  some  of  the 
ideas  underlying  reduced  integration. 

The  worked  example  considered  in  Section  2.3  used  a quadratic  shape  function.  The  equivalent  two-dimensional 
shape  function  would  be  a quadratic  Lagrange  shape  function  used  in  a rectangular  element.  The  use  of  an  isoparametric 
formulation  and  distorted  elements  probably  degrades  the  accuracy  to  be  expected  from  a four  point  Gauss  quadratiue 
formula.  However,  as  is  apparent  from  Section  3,  the  use  of  reduced  integration  produces  considerable  improvement. 
Similar  remarks  apply  to  the  Serendipity  quadratic  element. 

The  linear  rectangular  element  also  produces  more  accurate  results  when  reduced  integration  (one  sampling  point 
per  element)  is  used.  However  the  improvement  over  the  use  of  exact  integration  is  not  as  good  as  for  quadratic  rect- 
angular elements.  This  is  in  spite  of  the  fact  that  the  ratio  of  function  evaluation  points  per  unknown  is  the  same  as  for 
the  use  of  Lagrange  quadratic  rectangular  elements  (Tables  7 and  8).  Strang  and  Fix(ref.l3)  indicate  that  the  use  of  one 
integration  point  with  a linear  shape  function  over  a single  element  leads  to  an  indefinite  result,  and  caution  against  its 
use.  Even  though  typical  results  with  a linear  rectangular  element  and  reduced  integration  produced  a smaller  r.m.s. 
error  than  the  use  of  exact  integration  the  solution  was  often  physically  unsatisfactory  (e.g.  figure  14). 

The  use  of  reduced  integration  with  triangular  elements  appears  from  this  study  to  be  relatively  ineffective.  In  general 
terms,  since  triangular  elements  do  not  require  as  many  unknown  parameters  they  do  not  admit  as  many  higher  order  terms 
as  rectangular  elements,  and  hence  are  more  likely  to  be  integrated  with  low  order  formulae.  This  is  the  case  for  linear 
triangular  elements  where  even  one  sampling  point  per  element  is  sufficient  to  produce  exact  integration.  This  is 
supported  by  the  relatively  large  number  of  function  evaluation  points  per  unknown  (see  Tables  7 and  8).  The  case  of 
quadratic  triangular  elements  is  inconclusive.  The  use  of  a four  point  integration  formula  certainly  doesn’t  produce 
exact  integration.  However  the  interpretation  of  the  sampling  points  as  giving  rise  to  a special  fit  (e.g.  least-squares)  of 
the  residual  is  lacking  for  triangular  elements.  Consequently,  perhaps  it  is  not  surprising  that  the  use  of  a four  point 
integration  formula  has  not  produced  better  results  than  the  use  of  a seven  point  formula  on  average.  In  some  instances 
the  results  were  significantly  worse.  The  question  as  to  whether  a symmetric  four  point  formula  exists  that  possesses 
a comparable  interpretation  to  that  of  a four  point  Gauss  quadrature  formula  for  rectangular  elements,  deserves  further 
study. 

If  the  reason  for  the  success  of  reduced  integration  is  temporarily  disregarded  it  is  apparent  that  the  ability  of  rectangular 
elements  to  respond  to  the  use  of  reduced  integration  makes  them  significantly  more  efficient  than  triangular  elements. 

The  results  of  Section  3 indicate  that  the  quadratic  Serendipity  element  is  more  efficient  than  the  quadratic  Lagrange 
element  and  both  are  more  efficient  than  the  linear  element.  This  order  of  merit  is  consistent  with  the  results  of 
reference  5. 
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The  present  study  dues  not  alter  the  general  conclusion  of  reference  5 that  the  quadratic  triangular  element  is  more 
efficient  tlian  the  linear  triangular  clement. 


5.  CONCLUSIONS 

From  the  present  study  of  a Galerkin  finite  element  formulation  for  a representative  example  of  incompressible 
inviscid  flow,  the  following  conclusions  have  been  drawn: 

(a)  With  exact  numerical  integration  the  finite  element  solution  is  in  much  closer  agreement  with  the  exact 
solution  at  the  Gauss  integration  points  than  at  the  nodal  points. 

(b)  If  the  solution  at  the  Gauss  points  is  extrapolated  to  the  edge  of  the  element  the  solution  is  not  significantly 
more  accurate  than  the  solution  at  the  nodal  points. 

(c)  Results  for  a one-dimensional  problem  suggest  that  the  use  of  reduced  integration,  i.e.  a lower  order  Gauss 
quadrature  formula,  is  equivalent  to  fitting  the  equation  residual  in  a least-squares  sense  over  each  element. 

(d)  Reduced  integration  has  produced  results  that  are  as  accurate  as  those  produced  by  exact  integration 
requiring  ten  times  the  CPU  time,  when  used  with  quadratic  rectangular  elements  of  the  Serendipity  type. 

(e)  When  applied  to  quadratic  rectangular  elements  of  the  Lagrange  type,  the  use  of  reduced  integration  has 
produced  solutions  of  sufficient  accuracy  to  require  from  S to  10  times  as  much  CPU  time,  depending  on 
grid  refinement,  to  produce  comparable  results  when  exact  integration  is  used. 

(0  The  use  of  a four  point  integration  formula  has  produced  inferior  results  to  the  use  of  a seven  point  integration 
formula  when  quadratic  triangular  elements  are  used  on  a refined  grid.  However  the  converse  is  true  for  coarse 
grid. 

(g)  The  use  of  a one  point  integration  formula  with  linear  rectangular  elements  has  produced  superior  results  to 
the  use  of  a four  or  nine  point  integration  formula  with  the  same  elements. 

(h)  The  use  of  a one  point  integration  formula  with  linear  triangular  elements  has  produced  exactly  the  same 
result  as  the  use  of  a seven  point  integration  formula  with  the  same  elements. 

(i)  The  improvement  associated  with  the  use  of  reduced  integration  on  Serendipity,  rectangular  elements  has 
been  found  to  be  independent  of  grid  refinement. 

(j)  The  lack  of  success  of  reduced  integration  applied  to  quadratic  triangular  elements  is  probably  due  to  the 
failure  of  the  integration  formula  used  to  approximate  the  residual  in  the  least -squares  sense. 

(k)  The  effectiveness  of  quadratic,  rectangular.  Serendipity  elements  foi  the  use  of  reduced  integration  is  consistent 
with  comparable  results  obtained  from  structural  applications  of  the  finite  element  method. 

This  study  indicates  that  the  combination  of  reduced  integration  and  quadratic  Serendipity  elements  will  permit 
accurate  computational  results  to  be  obtained  for  flows  about  three-dimensional  bodies  with  much  coarser  grids,  and 
hence  smaller  execution  times,  than  would  have  been  possible  otherwise. 
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NOTATION 


A 

B. 


CPU 
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R 
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a.. 

‘1 


b.. 

•J 
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q 
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‘*T 

|ql 

u,  V 
U,  V 


area  of  triangle.  Appendix  1 


weighted  integral  over  known  q/s 


central  processing  unit 
term  in  the  stiffness  matrix 


number  of  nodal  unknowns  to  be  determined 
number  of  function  evaluations  per  equation 
total  number  of  function  evaluations  in  the  domain 


shape  function,  number  of  points  on  the  body 
residual 

weight  function 

coefficient  in  the  solution  matrix 

coefficient  in  the  governing  algebraic  equation,  Section  3 

number  of  contributing  elements 
total  number  of  elements  in  the  domain 

number  of  function  evaluations  per  element;  order  of  differentiation 
number  of  points  in  the  integration  formula;  number  of  nodal  unknowns 
number  of  nodes 


order  of  finite  element  representation 

exact  tangential  velocity  component  at  the  body  surface 

typical  nodal  unknown 

tangential  velocity  component  at  the  body  surface  calculated  from  the  finite  element  solution 

vector  of  nodal  unknowns 
velocity  components,  figure  1 
nodal  values  of  velocity  components 

weight  associated  with  the  function  value  at  the  k'^  point  in  the  numerical  integration  formula 
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x.y 


Vi-f, 


cartesian  coordinates 

root  mean  square  dilTerence  between  and  q^ 

angular  position  around  a circular  cylinder,  measured  from  the  front  stagnation  point 
triangular  coordinates,  Appendix  I 
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APPENDIX  I 

NUMERICAL  INTEGRATION  FORMULAE  FOR  VARIOUS  SHATC 
FUNCTION/ELEMENT  COMBINATIONS 


The  numerical  integration  over  a triangle  of  area  A is  obtained  from 


(II) 


where  and  f . are  triangular  coordinates  associated  with  the  i^^  sampling  point  and  w.  is  the  associated  weight. 

The  coordinates  and  weights  of  the  various  triangular  integration  formulae  are  set  our  below.  The  various  sampling 
points  are  shown  in  figure  20 
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y 


one  point, 

a = {1/3,  1/3, 1/3  1, 

W * 

a 

1 

four-point, 

a = {1/3, 1/3, 1/3  j. 

w * 
a 

-27  / 48 

b = {3/5, 1/5, 1/5  j. 

25/48 

c = {l/5,3/5,l/5  }, 

w * 
c 

25/48 

d ={  1/5,  1/5,  3/5  j. 

'"d  = 

25/48 

seven-point. 

a = {1/3, 1/3, 1/3  j. 

w = 
a 

0.225 

b = {ai,/3,,ft  1 , 

0.13239 

41528 

c = {|3,,a,,0,  j , 

w = 
c 

0.13239 

41528 

d ={  A . /3i . a,  j , 

'"d  = 

0.13239 

41528 

* “ {tt2, A,  A 1 , 

w » 
e 

0.12593 

91805 

I ~ j > 

Wf  = 

0.12593 

91805 

g = {A.A.aai  . 

w = 
g 

0.12593 

91805 

where 


ai 

= 0.05971 

58718 

A 

= 0.47014 

20641 

<h 

= 0.79742 

69854 

A 

= 0.10128 

65073 

The  seven-point  formula  has  been  used  to  produce  exact  integration  for  the  linear  and  quadratic  triangular  isoparametric 
elements.  Tire  four-point  formula  has  been  used  to  cany  out  reduced  integration  of  quadratic  triangular  isoparametric 
elements.  The  one-point  formula  has  been  used  to  cany  out  reduced  integration  of  linear  triangular  elements.  However 
the  one-point  formula  produces  exact  integration  for  this  problem. 
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Gauss  quadrature  formulae  of  the  form 
1 1 


J I d{.d,?=  ^ ^H..H..f(^.»j.) 


(1.2) 


i=l  j=l 


have  been  used  to  integrate  f over  rectangular  elements,  i , rj.  are  the  coordinates  of  the  i,  i‘"  point  and  H , H 

1 J 1 j 

are  the  corresponding  weights.  The  coordinates  and  weights  for  Gauss  quadrature  formulae  up  to  n = 3 are  set  out 
below.  The  various  sampling  points  are  shown  in  Hgure  20. 


one-point  (n=l) 

= 0.0 

H.,H.  = 2.0 
1 J 

four-point  (n=2) 

H..H.  = 1.0 
' J 

nine-point  (n=3) 

= ±V3/V^ 

H.,H.  = 5/9 
• J 

= 0.0 

= 8/9 

The  nine-point  formula  has  been  used  to  produce  exact  integration  for  quadratic  rectangular  isoparametric  elements 
of  both  the  Serendipity  and  Lagrange  kind.  The  four -point  formula  has  been  used  to  carry  out  reduced  integration 
over  quadratic  rectangular  isoparametric  elements  and  exact  integration  over  linear  rectangular  elements.  The  one 
point  formula  has  been  used  to  carry  out  reduced  integration  over  linear  rectangular  elements. 

The  definition  of  the  coordinates  used  in  the  triangular  and  rectangular  elements  and  a fuller  description  of  numerical 
integration  may  be  found  in  reference  1 1 . 
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TAbLL  3.  EXACT  INTEGRATION  FORMULAE  ON  COARSE  GRIDS 


Element 

type 

Shape 

function 

Number  of 
unknowns 

CPU 
time  (s) 

Nodal  r.m.s. 
difference 

Figure 

reference 

Rectangular 

Linear 

199 

5.1 

0.043 

Rectangular 

Quadratic  (S) 

149 

4.5 

0.049 

6 

Triangular 

Linear 

279 

8.5 

0.040 

Triangular 

Quadratic 

199 

7.1 

0.120 

Triangular 

Linear 

199 

5.1 

0.061 

Rectangular 

Quadratic  (L) 

199 

7.5 

0.047 

Rectangular 

Quadratic  (L) 

1S9 



5.0 

0.071 

S = Serendipity  L = Lagrange 


TABLi;  4.  REDUCED  INTEGRATION  FORMULAE  ON  COARSE  GRIDS 


Element  type  ! 
type 

Shape 

function 

Number  of 
unknowns 

CPU 
time  (s) 

Nodal  r.m.s. 
difference 

Figure 

reference 

Rectangular 

Linear 

199 

5.3 

0.033 

9.11 

Rectangular 

Quadratic  (S) 

149 

4.5 

0.015 

6,8,9,13,19 

Triangular 

Linear 

279 

8.5 

0.040 

10 

Triangular 

Quadratic 

199 

7.3 

0.060 

10,12 

Triangular 

Linear 

199 

5.1 

0.061 

11 

Rectangular 

Quadratic  (L) 

199 

7.7 

0.022 

12 

Rectangular 

Quadratic  (L) 

159 

5.1 

0.042 

13 
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TABLE  5.  EXACT  INTEGRATION  FORMULAE  ON  MODERATE  GRIDS 


Element 

type 

Sliape 

function 

Number  of 
unknowns 

1 

CPU 
time  (s) 

Nodal  r.m.s. 
difference 

Figure 

reference 

Rectangular 

Linear 

399 

18 

0.034 

Rectangular 

Quadratic  (S) 

299 

15 

0.023 

7 

Triangular 

Linear 

599 

34 

0.024 

Triangular 

Quadratic 

399 

23 

0.018 

Triangular 

Linear 

399 

17 

0.034 

Rectangular 

Quadratic  (L) 

399 

25 

0.023 

Rectangular 

Quadratic  (L) 

319 

16 

0.025 

Rectangular 

Quadratic  (S) 

582 

65 

0.013 

TABLE  6.  REDUCED  INTEGRATION  FORMULAE  ON  MODERATE  GRIDS 


Element 

type 

Shape 

function 

Number  of 
unknowns 

CPU 
time  (s) 

Nodal  r.m.s. 
difference 

Figure 

reference 

Rectangular 

Linear 

399 

17 

0.022 

14, 16 

Rectangular 

Quadratic  (S) 

299 

15 

0.0071 

7,14,18,19 

Triangular 

Linear 

599 

33 

0.024 

15 

Triangular 

Quadrat  ic 

399 

22 

0.026 

15,17 

Triangular 

Linear 

399 

16 

0.034 

16 

Rectangular 

Quadratic  (L) 

399 

25 

0.0072 

17 

Rectangular 

Quadratic  (L) 

319 

17 

0.0080 

18 

Rectangular 

Quadratic  (S) 

582 

69 

0.0035 

19 
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TABLE  7.  DEGREE  OF  DEFINITENESS  FOR  COARSE  MESH 


Reduced  integration 

Exact  integration 

Element/ 
shape  function 
combination 

Number  of 
unknowns 

L 

Number  of 
elements 

e 

No.  of 
function 
eval. 
points 

Function 
eval.  per 
unknown 

Hp/L 

Number  of 
unknowns 

L 

Number  of 
elements 

e 

No.  of 
function 
eval.. 
points 

Mt 

Function 
eval.  per 
unknown 

Mj./L 

Linear  rectangle 

199 

100 

200 

1.005 

199 

100 

800 

4.020 

Quad,  rectangle 
(Serendipity) 

149 

25 

200 

1.342 

149 

25 

450 

3.020 

Quad,  rectangle 
(Lagrange) 

199 

25 

200 

1.005 

199 

25 

450 

2.261 

Linear  triangle 

199 

200 

400 

2.010 

199 

200 

1600 

8.040 

Quad,  triangle 

199 

50 

400 

2.010 

199 

50 

700 



3.518 

TABLE  8.  DEGREE  OF  DEFINITENESS  FOR  MODERATE  MESH 


Reduced  integration 

Exact  integration 

Element/ 
shape  function 
combination 

Number  of 
unknowns 

L 

Number  of 
elements 

e 

No.  of  j 
function 
eval. 
points 

■^T 

Function 
eval.  per 
unknown 

M.j,/L 

Number  of 
unknowns 

L 

Number  of 
elements 

e 

No.  of 
function 
eval. 
points 

M.J. 

Function 
eval.  per 
unknown 

M.J./L 

Linear  rectangle 

399 

200 

1 

400 

1.003 

399 

200 

1600 

4.010 

Quad,  rectangle 
(Serendipity) 

299 

50 

400 

1.338 

299 

50 

900 

3.010 

Quad,  rectangle 
(Lagrange) 

399 

50 

400 

1.003 

399 

50 

900 

2.256 

Linear  triangle 

399 

400 

800 

2.005 

399 

400 

3200 

8.020 

Quad,  triangle 

399 

100 

800 

2.005 

399 

100 

1400 

3.509 
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Quadratic  rectangular  elements  (Serendipity) 

Solution  at  integ.  pts.  extrapolated  to 
body  surface  (Solution  C,  Table  2) 

O Nodal  pt.  solution  (Solution  B,  Table  2) 

^ Nodal  pt.  solution  extrapolated  to  the  same  points  as 
extrapolated  integ.  pts.  solution  (Solution  D,  Table  2) 


Tangential  velocity  difference  (q 


WRE-TR-I810(W) 
Figure  S 


Quadratic  rectangular  elements  (Serendipity) 

Q Solution  at  integ.  pts.  extrapolated  to 
body  surface  (Solution  C,  Table  2) 

O Nodal  pt.  solution  (Solution  B,  Table  2) 

^ Nodal  pt.  solution  extrapolated  to  the  same  points  as 
extrapolated  integ.  pts.  solution  (Solution  D,  Table  2) 


Figure  S.  Variation  of  tangential  velocity  differences  with  angular  position  : 9 point  integration  formula 


Tangential  velocity. 
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CPU  = 15s 

O 4 pt.  integ.formula  | quadratic  rectangular 
13  9 pt.  integ.formula  [ element  (Serendipity) 
— — — Exact 

r.m.s.  error  (4  pt.  i.f.)  = 0.(X)7 
r.m.s.  error  (9  pt.  i.f.l  = 0.023 


0 10  20  30 

Angular  position,  6 (degrees) 

Figure  7.  Comparison  of  accuracy  of  reduced  integration  and  exact  integration  on  a moderate  grid 


Tangential  velocity. 


Figure  9.  Comparison  of  linear  and  quadratic  shape  functions  - rectangular  elements  ■ coarse  grid  - reduced  integration 
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Tangentical  velocity. 
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A Linear  triangular  element 

0 Linear  rectangular  element 

— — Exact  solution 


Angular  position,  6 (degreees) 


Figure  1 1 . Comparison  of  triangular  and  rectangular  elements  - linear  shape  function  - coarse  grid  • reduced  integration 
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A Linear  triangular  element 

A Quadratic  triangular  element 


Angular  position,  0 (degrees) 


Figure  I 5.  Compurison  of  linear  and  quadratic  shape  functions  - triangular  elements  - moderate  grid  • reduced  integration 


X-  Function  evaluation  points 
for  numerical  integration 


Tangential  velocity. 
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One  - Point  Four  ■ Point  Seven  - Point 

(a)  Sample  points  for  triangular  elements 


One  - Point  Four  ■ Point  Nine  - Point 

(b)  Sample  points  for  rectangular  elements 


Figure  20.  Sample  points  for  numerical  integration  formulae 
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